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Hydro  instabilities  have  been  identified  as  a  potential  cause  of  performance  degradation  in  inertial 
confinement  fusion  (ICF)  experiments.  We  study  instabilities  associated  with  a  single  Richtmyer 
—Meshkov  (RM)  interface  in  a  circular  geometry,  idealized  from  an  ICF  geometry.  In  an  ICF  applica¬ 
tion,  atomic  level  mix,  as  an  input  to  nuclear  burn,  is  an  important,  but  difficult  to  compute,  variable.  We 
find  numerical  convergence  for  this  important  quantity,  in  a  purely  hydro  study,  with  only  a  mild 
dependence  on  the  Reynolds  number  of  the  flow,  in  the  high  Reynolds  number  limit.  We  also  find  that 
mixing  properties  show  a  strong  sensitivity  to  turbulent  transport  parameters;  this  sensitivity  translates 
into  an  algorithmic  dependence  and  a  nonuniqueness  of  solutions  for  nominally  converged  solutions.  It  is 
thus  a  complication  to  any  verification  and  validation  program.  To  resolve  the  nonuniqueness  of  the 
solution,  we  propose  a  validation  program  with  an  extrapolation  component,  linking  turbulent  transport 
quantities  in  experimental  regimes  to  mildly  perturbed  turbulent  transport  values  in  ICF  Reynolds 
number  regimes.  In  view  of  the  observed  solution  nonuniqueness,  the  validation  program  and  its  jus¬ 
tification  from  the  results  presented  here,  has  a  fundamental  significance. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

1.1.  Background 

Hydro  instabilities  are  known  to  limit  inertial  confinement  fusion 
(ICF)  neutron  production  [  1  ].  See  also  the  more  recent  discussions  in 
Ref.  [2].  For  predictions  of  neutron  production,  the  atomic  level 
description  of  temperature  and  mixture  concentrations  provide  the 
main  input  to  a  thermonuclear  (TN)  burn  computation  [3].  The  in¬ 
terest  in  mix  is  brought  about  by  current  experiments  of  ICF  pro¬ 
cesses.  For  instance,  there  is  an  observed  degradation  of  40-75%  in 
neutron  production,  attributed  to  atomic  mixing,  estimated  [4] 
based  on  experiments  conducted  on  the  Omega  laser.  In  addition, 
it  is  proposed  that  low  mixing  is  required  to  achieve  ignition,  spe¬ 
cifically,  the  required  maximum  level  of  contaminated  fuel  from  the 
ablator  material  must  be  less  than  an  estimated  25%-40%  [1  ].  Modal 
growth  factors  of  10—100  or  more  are  reported  in  the  linearized 
analysis  of  [  1  ],  suggesting  a  role  for  the  nonlinear  mode  coupling  and 
even  perhaps  the  fully  developed  turbulent  mixing  regimes,  beyond 
the  ablation  stabilized  weakly  nonlinear  theory  reviewed  in  Ref.  [5]. 
ICF  simulations  do  not  predict  correctly  the  observed  dynamics  of 
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ICF  capsules,  which  exhibit  significant  yield  degradation  relative  to 
design.  Thus  a  re-examination  of  the  validation  process  would 
appear  to  be  appropriate.  Experimental  analysis  [6]  of  the  ICF  hot 
spot  mix  is  based  on  Ge  doped  ablator  material  found  within  the  hot 
spot  at  ignition  time.  Only  mix  at  a  level  as  predicted  by  design 
analysis  was  observed. 

Known  ICF  instabilities  or  asymmetries  occur  primarily  in  four 
distinct  stages  during  an  ICF  process.  They  begin  with  (1)  the  laser 
drive  and  are  followed  by  (2)  instability  growth  at  the  ablation 
surface.  In  addition,  (3)  Richtmyer— Meshkov  (RM)  instabilities 
occur  at  the  boundary  of  ablator  to  deuterium-tritium  (DT)  ice  and 
the  boundary  of  DT  ice  to  DT  gas.  These  propagate  inward  to  define 
initial  conditions  for  the  (4)  deceleration  phase,  which  is  known  to 
be  strongly  RT  unstable.  The  second  of  these  instabilities  is  an 
ablation  modified  RT  instability,  the  third  is  RM  and  the  fourth  is  RT 
in  nature.  In  Ref.  [7],  velocity  perturbations  of  1%  at  the  beginning  of 
the  deceleration  phase  (4)  are  considered  to  be  dangerous  for  the 
performance  of  the  capsule.  The  Verification  and  Validation  (V&V) 
methods  of  this  paper  are  applicable  to  the  second,  third  and  fourth 
instabilities  mentioned  above.  Inaccuracies  in  prediction  of  the 
yield  cliff  suggests  a  possible  need  for  improved  modeling  of  the 
ablation  phase,  such  as  the  inclusion  of  subgrid  scale  (SGS)  models 
as  proposed  in  Ref.  [2],  and  the  value  in  reexamination  of  any  of  the 
four  instability  stages,  as  well  as  other  issues  not  discussed  here. 
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Here  we  address  V&V  aspects  of  a  hydro  instability  study  related 
to  these  broad  ICF  concerns.  We  study  the  RM  instability,  which  is 
a  primary  ICF  hydro  instability,  in  its  relation  to  atomic  scale  mix; 
the  conclusions  have  implications  for  RT  V&V  as  well. 

The  main  thrust  of  the  paper  is  to  offer  scientific  data  in  support 
of  a  proposal  for  validation  with  extrapolation.  We  use  a  single  code 
with  no  free  parameters,  validated  in  common  for  RT  and  RM  data. 
This  goal,  infeasible  for  Reynolds  Averaged  Navier- Stokes  (RANS),  is 
a  natural  objective  for  Large  Eddy  Simulations  (LES),  which  use  dy¬ 
namic  SGS  models.  Validation  for  RT  was  previously  considered  and 
is  summarized  in  Sec.  1.4.  The  RT  case  is  more  sensitive  to  transport 
(molecular,  turbulent  or  numerical)  than  are  RM  instabilities.  RT  is 
thus  more  demanding  for  validation.  This  difficulty  can  be  inferred 
from  the  fact  that  (in  the  authors’  opinion)  the  FT/LES/SGS  code 
discussed  here  appears  to  be  unique  as  a  multimode  RT  validated 
compressible  LES  code.  For  RM,  we  rely  on  an  extensive  code  com¬ 
parison  study  [8]  with  RAGE  [9].  Generally  excellent  agreement  was 
found  for  most  measures  of  comparison,  however,  exceptions  were 
noted  for  the  temperature  field  (resolved  in  favor  of  our  code  Fron- 
Tier  (FT)).  In  the  case  of  high  density  ratios,  such  as  those  that  occur 
at  the  hot/cold  DT  interface  in  an  ICF  capsule,  the  edges  of  the  mixing 
zone  showed  slow  numerical  convergence,  especially  for  RAGE, 
indicating  potential  problems  for  commonly  used  resolution  levels. 
Translated  into  the  ICE  context,  this  raises  the  question  of  simulation 
velocity  perturbations  at  the  beginning  of  the  deceleration  phase 
and  perhaps  of  penetration  of  cold  DT  gas  from  the  original  DT  ice 
layer  of  an  ICF  capsule  into  the  ICE  hot  spot,  or  of  ablator  material  (as 
has  been  observed  experimentally),  any  of  which  can  potentially 
degrade  performance.  Here  we  are  mainly  concerned  with  the 
extrapolation  from  the  experimental  regime  to  an  ICE  regime  or  to 
still  higher  Reynolds  numbers.  For  the  purpose  of  this  paper,  we 
understand  the  experimental  regime  to  be  Re  ~  3.5  x  10^,  typical  of 
laboratory  RT  experiments  and  the  ICF  regime  to  be  Re  ~  6  x  10^, 
with  higher  values  of  Re  also  considered.  We  observe  that  the 
obvious  step  of  RT  validation  using  NIF  or  Omega  laser  data  does  not 
address  the  multimode,  mode  coupling  RT  growth  stage,  as  the  laser 
experiments  performed  to  date  have  not  reached  this  time  regime. 

The  two  main  steps  suggested  here  (RT  validation  and  assess¬ 
ment  of  high  Re  extrapolation,  with  a  zero  parameter,  theoretically 
constrained  model)  are  not  part  of  common  V&V  methodology. 
Most  multiphysics  codes  fail  for  RT  validation  in  the  LES  regime, 
and  generally  ICE  codes  do  not  employ  dynamic  subgrid  scale 
models,  which  allow  for  smooth  extrapolation  from  experimental 
to  ICF  Reynolds  numbers.  The  validation  proposed  here  is  suggested 
as  a  complement  to,  or  preliminary  to,  and  not  a  replacement  for, 
validation  as  presently  performed.  The  most  fundamental  vali¬ 
dation  (as  yet  apparently  still  an  open  question)  is  comparison  of 
simulations  against  ICE  capsule  performance  directly. 

The  issue  we  raise  is  that  the  atomic  mix  properties  are  sensitive 
to  the  turbulent  transport  coefficients,  and  through  them  to  the 
numerical  algorithm.  Stated  more  directly,  apparently  converged 
solutions  are  nonunique.  They  depend  in  an  essential  manner  on 
details  of  the  numerical  algorithm.  This  unfortunate  fact  is  docu¬ 
mented  here.  However,  we  find  only  a  mild  variation  in  these  co¬ 
efficients  in  passing  from  an  experimental  regime  to  the  ICF  regime, 
as  far  as  purely  hydro  issues  are  concerned.  Thus  a  mild  extrap¬ 
olation  from  experimentally  validated  transport  coefficients  is 
feasible,  yielding  a  significant  constraint  on  the  turbulent  transport 
terms.  We  bridge  this  gap  using  algorithmic  choices  supported  and 
constrained  by  theory  and  mathematical  verification.  We  propose 
a  program  to  carry  out  this  objective,  based  on  (a)  a  numerical 
algorithm  with  front  tracking,  to  limit  numerical  thermal  conduc¬ 
tion  and  concentration  diffusion  and  on  (b)  dynamic  subgrid 
models  to  assure  a  theoretically  consistent  parameter  free  setting  of 
the  turbulent  transport  coefficients. 


The  context  of  the  present  study  is  a  circular  Richtmyer— 
Meshkov  fluid  instability  problem;  in  view  of  the  twin  computa¬ 
tional  requirements  of  extreme  levels  of  mesh  refinement  and  of 
a  parametric  variation  study  involving  multiple  simulations,  we 
consider  this  problem  in  2D.  At  high  Re,  with  ionized  molecules,  we 
consider  a  fluid  with  plasma  like  properties.  These  particles  may  be 
strongly  coupled  to  a  radiation  held,  while  the  thermal  conduction 
between  ions  may  be  mediated  by  the  ion-electron  thermal  cou¬ 
pling  and  the  faster  electron  thermal  diffusion.  To  model  at  a  purely 
hydro  level,  we  assume  a  highly  heat  conductive  Prandtl  number 
Pr  =  10“"^  and  a  Schmidt  number  Sc  =  1.  We  define  the  high  Re  limit 
to  be  taken  as  viscosity  v  ^  0,  with  fixed  values  for  Sc  and  Pr.  The 
parameter  study  analyzes  a  series  of  mesh  levels  and  Re  values,  up 
to  Re  =  6  X  10^  and  Re  =  00  ,  extending  a  study  [10]  of  the  same 
problem.  In  future  work,  we  will  include  plasma  and  radiation 
effects. 

1.2.  Summary  of  main  results 

Here,  we  explain  a  fundamental  V&V  challenge  and  present 
scientific  data  which  supports  a  V&V  route  to  bypass  these  obsta¬ 
cles.  The  challenge  is  the  nonuniqueness  and  algorithmic  depend¬ 
ence  of  nominally  converged  solutions  for  turbulent  mixing 
problems.  The  route  to  bypass  this  obstacle  is  the  combined  use  of 
parameter  free  dynamic  SGS  models  and  control  over  numerical 
diffusion  through  front  tracking. 

To  validate  this  algorithmic  choice,  we  use  a  two  step  V&V 
method.  The  main  step  is  RT  validation  against  the  wealth  of  RT 
experimental  data  at  Re  ~  3.5  x  10^.  This  data  is  highly  sensitive  to 
transport  (molecular,  turbulent,  numerical)  to  the  point  that  the  FT/ 
LES/SGS  code  used  here  appears  to  be  the  unique  compressible 
multiphysics  LES  code  to  have  been  validated  in  a  zero  parameter 
manner  against  this  data.  This  first  step  was  reported  in  earlier 
publications  [11,12]  and  is  summarized  here. 

The  second  step  in  our  V&V  method  is  an  extrapolation  from 
Re  ~  3.5  X  10"^  to  the  ICE  regime  of  Re  ~  6  x  10^  or  higher  up  to 
Re  =  00 .  This  step  leads  to  the  main  technical  result  of  this  paper.  In 
essence,  we  show  (a)  that  the  consequences  of  the  extrapolation 
are  not  large  in  their  influence  on  numerical  simulations,  and  (b) 
that  there  is  an  available  theory  to  constrain  the  extrapolation.  The 
theory  [13],  which  is  the  basis  of  dynamic  SGS  models,  is  not  new  to 
this  paper. 

Thus  our  main  results  can  be  stated  as  (a)  numerical  con¬ 
vergence  of  the  turbulent  mixing  statistics  in  question  and  (b) 
a  small  magnitude  of  the  perturbation  in  these  statistics  due  to 
a  change  of  Re. 

1.3.  The  numerical  algorithm 

The  numerical  study  is  based  on  the  front  tracking  code  FronTier 
[14,15],  enhanced  with  an  LES  turbulence  model  having  dynamic 
SGS  terms  [13].  We  refer  to  this  algorithm  as  FT/LES/SGS.  These  two 
algorithms  are  built  upon  a  hydro  package  (compressible,  with 
general  equation  of  state)  using  the  MUSCL  [16]  algorithm.  To 
facilitate  coupling  of  these  features  with  other  hydro  packages,  we 
have  developed  an  API.  To  avoid  complications  with  multispecies 
diffusion  [17],  we  assume  mixing  of  two  fluids  only. 

1.3.1.  Front  tracking 

With  front  tracking,  our  aim  is  to  achieve  enhanced  resolution  of 
interfaces  between  fluids  and  maintain  this  resolution  throughout 
the  simulation.  This  allows  us  to  minimize  the  numerical  mass 
diffusion  associated  with  steep  concentration  or  thermal  gradients. 
We  see  improved  results  even  for  coarser  grids,  since  the  dis¬ 
continuity  that  normally  leads  to  the  smearing  in  coarser  grids  is 
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dynamically  tracked,  preventing  one  fluid  from  numerically  inter¬ 
acting  across  the  interface  with  the  other. 

We  accomplish  this  goal  by  creating  a  separate,  lower  dimen¬ 
sional  mesh  along  the  front  (interface)  which  separates  the  fluids, 
for  example  following  an  iso-temperature  contour  within  a  sharp 
thermal  gradient.  The  front  then  separates  fluids  with  distinct 
physical  characteristics.  During  each  time  step  we  propagate  the 
front’s  mesh  to  its  new  location  based  on  the  surrounding  interior 
states.  The  front  points  are  propagated  in  a  normal  direction,  uti¬ 
lizing  a  Riemann  problem  formulation  to  resolve  the  different  ve¬ 
locities  obtained  via  one  sided  extrapolation  from  interior  states  to 
the  front  location. 

When  interfacial  tangling  occurs,  we  utilize  algorithms  to  un¬ 
tangle  and  redistribute  points  on  the  front  mesh,  allowing  for 
a  logical  reconstruction  of  a  tangled  interface.  For  more  information 
on  this  process  and  the  FronTier  front  tracking  algorithm,  we  refer 
the  reader  to  [15,18,19]. 

The  front  mesh  is  overlayed  on  a  normal  Eulerian  grid.  Where 
it  intersects  cells,  we  cut  those  cells  into  multiple  pieces  to 
ensure  separation  of  the  distinct  fluids  across  an  interface, 
illustrated  in  Fig.  1. 

When  updating  interior  states  on  a  cut  cell,  ghost  cells  con¬ 
structed  from  states  defined  in  neighboring  cells  on  the  same  side 
of  the  front  allow  for  a  1 -sided  update,  consisting  only  of  states 
from  those  cells  on  the  same  side  of  the  discontinuity  interface. 

This  construction  allows  us  to  control  numerical  diffusion, 
especially  with  complicated  interfaces,  as  arise  in  hydro  in¬ 
stabilities  in  ICF  processes. 

13.2.  LES  with  dynamic  subgrid  models 

The  purpose  of  subgrid  models  is  to  capture  the  effect  of  the 
unresolved  scales  (below  the  grid  level)  on  those  that  are  resolved. 
We  start  with  the  interpretation  that  the  numerical  solution  values 
represent  grid  cell  averaged  quantities.  Conventionally,  in  the  der¬ 
ivation  of  Altered  Navier-Stokes  equations  for  use  in  Large  Eddy 
Simulations,  there  is  a  convolution  average  by  a  positive  weight 
function  (the  Alter),  and  the  equations  are  derived  for  the  Altered 
quantities.  We  omit  this  Alter  step,  and  use  the  cell  averaged 
quantities  directly.  The  average  of  the  nonlinear  terms  introduces 
new  unknowns  into  the  equations  for  which  new  equations,  called 
closures,  are  required.  These  show  up  as  flux  type  terms,  and  are 
the  origin  of  the  turbulent  transport.  The  closure  term  must  be 
modeled;  often  a  solution  gradient  is  used  for  a  turbulent  flux,  with 
an  undetermined  coefficient.  With  dynamic  subgrid  models 
[13,20,21],  the  otherwise  missing  coefficient  is  determined  by  the 
numerical  solution  itself,  and  varies  with  space-time  according  to 
local  flow  conditions.  The  determination  is  achieved  by  consid¬ 
eration  of  the  closure  on  the  current  grid  and  on  a  once  coarser  grid. 


and  by  an  asymptotic  assumption  that  the  model  coefficient  is 
a  mesh  convergent  or  asymptotically  mesh  independent  quantity, 
so  that  a  common  value  can  be  used  between  the  two  grid  levels 
considered.  The  advantage  of  cell  averages  as  opposed  to  a  Alter  is 
the  absence  of  loss  of  spatial  resolution  and  of  solution  fluctuations 
averaged  over  by  the  Alter. 

1.3.3.  API 

The  API  is  presently  under  development  http://www.ams. 
sunysb.edu/  --  rkaufman/api. 

Presently,  it  links  a  client  code  to  the  FT/LES/SGS  front  tracking 
library.  All  functions  in  the  API  are  designated  as  belonging  to  the 
API  or  to  the  client.  Those  belonging  to  the  client  depend  on  client 
specific  information,  such  as  the  data  layout  of  solution  state  values. 
For  these  functions,  the  API  includes  a  reference  implementation  in 
a  few  regular  cases,  such  as  for  regular  rectangular  grid  based  state 
data.  The  main  point  of  the  API  functions  is  to  communicate  data 
between  the  client  and  the  front  library.  This  communication  is,  of 
course,  two  directional.  From  the  client  data,  we  extract  velocities, 
used  to  advance  the  front  in  time.  From  the  front  library,  we  extract 
front  location  information,  used  in  a  client  function  to  label  any 
given  flnite  difference  stencil  as  regular  or  irregular,  according  to 
whether  the  stencil  does  not  or  does  cross  a  tracked  front.  Finally, 
a  client  function  will  offer  modifled  (tracked)  stencil  state  values,  to 
affect  ghost  cells  in  a  neighborhood  of  the  front.  The  modifled 
stencil  of  state  values  will  be  differenced  using  the  usual  client 
difference  algorithm. 

1.4.  Prior  validation  and  code  comparison  studies  for  FT/LES/SGS 

1.4.1.  Rayleigh-Taylor 

FT/LES/SGS  has  achieved  systematic  agreement  with  a  wide 
range  of  Rayleigh-Taylor  experiments,  including  mixing  of  both 
concentration  and  thermal  flelds.  Uncertainty  over  initial  condi¬ 
tions  (not  recorded)  has  been  resolved,  by  backward  solution  of  the 
fluid  equations  from  an  early  time  (recorded)  to  an  initial  time  (not 
recorded),  and  with  uncertainty  quantiflcation  for  the  possible  er¬ 
rors  in  the  reconstruction  [22].  Agreement  is  precise  enough  to 
distinguish  between  distinct  experiments  and  their  distinct  values 
of  a.  See  Table  1. 

1.4.2.  Richtmyer-Meshkov 

FronTier  was  used  in  a  three  way  code  comparison  (with  RAGE 
and  PROMETFIEUS),  with  experiment  and  with  theory  [26]  for 
single  mode  Richtmyer-Meshkov  instabilities.  Extensive  compar¬ 
ison  of  FT/LES/SGS  with  the  already  validated  RAGE  code  for  mul¬ 
timode  Richtmyer-Meshkov  instabilities,  before  and  after  reshock, 
gave  generally  excellent  results  [8],  with  exceptions  related  to  the 


Fig.  1.  Discretized  interface  and  cut  cell  construction.  Left:  Tracked  front  overlayed  over  Eulerian  grid,  discretized  at  front  points.  Right:  Display  of  cut  cells  created  when  an  interface 
cuts  through  a  cell. 
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Table  1 

Experimental  and  FT/LES/SGS  simulation  values  of  a  compared  for  a  series  of  ex¬ 
periments,  with  error  bars,  if  provided,  by  either  the  simulation  or  the  experiment, 
and  the  discrepancy  if  any. 


Experiment 

Q!exp 

^^sim 

Discrepancy 

Smeeton-Youngs  [23]  #112 

0.052 

0.055 

6% 

Smeeton-Youngs  [23]  #105 

0.072 

0.072-0.080 

0% 

Smeeton-Youngs  10  exp. 

0.055-0.077 

0.066 

0% 

Andrews— Banerjee  [24] 

0.065-0.07 

0.069 

0% 

Andrews— Muschke  [25] 

0.059-0.081 

0.075 

0% 

Andrews-Muschke  [25] 

0.08-0.09 

0.084 

0% 

thermal  field  as  noted  above.  In  view  of  the  more  difficult  RT 
thermal  validation  and  the  absence  of  thermal  RM  experimental 
data,  we  plan  to  proceed  to  the  next  stage  of  our  studies,  which  is  to 
bring  the  simulations  closer  to  an  ICF  experimental  plane. 

1.5.  Outline  of  paper 

In  Sec.  2,  we  expand  on  a  method  for  stochastic  convergence  of 
turbulent  solutions  we  previously  introduced  and  lay  the  frame¬ 
work  for  the  convergence  strategy  of  our  algorithm.  In  Sec.  3,  we 
study  the  convergence  of  the  solution  PDFs  and  their  variation  with 
Re.  In  Sec.  4,  we  show  that  the  PDFs  are  sensitive  to  the  definition  of 
turbulent  transport.  In  other  words,  we  establish  numerical  non¬ 
uniqueness  for  turbulent  mixing  in  the  high  Re  limit.  In  Sec.  5  we 
outline  our  ideas  for  the  validation  of  numerical  studies  of  hydro 
instabilities  related  to  ICF  performance. 

2.  Stochastic  convergence  of  turbulent  solutions 

2.1.  Convergence  analysis  method 

Our  convergence  strategy  is  stochastic  in  nature  and  it  is  algo¬ 
rithmically  a  post  processing  method.  Therefore,  it  does  not  modify 
the  simulation  itself.  The  authors  are  not  aware  of  a  similar  con¬ 
struction  employed  by  others  for  analysis  of  numerical  solutions  of 
partial  differential  equations.  We  retain  solution  stochastic  infor¬ 
mation  (fluctuations  and  their  PDFs)  through  a  coarse  graining 
strategy  where  we  organize  the  mesh  cells  into  a  coarse  grained 
mesh  of  supercells.  All  solution  values  in  a  single  supercell  are 
combined  to  form  a  flnite  approximation  to  a  probability  density 
function  (PDF)  or  its  indeflnite  integral,  the  cumulative  distribution 
function  (CDF).  In  this  construction,  the  coarse  grained  supercell 
provides  the  (reduced)  spatial  or  space-time  localization.  Such 
a  construction  is  a  discrete  approximation  to  a  space-time 
dependent  probability  measure,  also  known  as  a  Young  Measure. 


Whether  assessing  convergence  or  the  effects  of  variation  of 
physical  parameters,  we  are  comparing  solutions.  The  difference 
between  solutions  is  expressed,  in  each  supercell,  as  the  difference 
between  two  coarse  grain  localized  CDFs.  In  each  supercell,  the 
solution  data  is  discretized  into  bins  based  on  the  solution  variables 
under  study.  From  this  binned  data,  we  create  an  empirical  PDF  or 
CDF.  The  norm  of  the  difference  is  deflned  first  through  an  Li  norm 
over  the  solution  state  variables,  i.e.,  of  the  Li  norm  difference  of  the 
CDFs  over  the  space  of  bins,  and  normalized  by  division  by  the 
number  of  bins,  and  by  division  by  the  difference  between  the 
maximum  and  minimum  solution  values  over  the  domain.  Such 
a  norm  difference  is  deflned  as  a  function  deflned  on  each  supercell. 
A  spatial  mean  (an  Li  norm  over  all  supercells  in  the  domain, 
divided  by  the  total  supercell  area  in  the  domain)  completes  the 
deflnition  of  the  norm  of  the  difference  of  the  CDFs.  This  strategy 
for  convergence  has  been  developed  in  previous  papers  [11,27,28]. 

2.2.  Supercell  and  bin  sizing 

We  now  examine  choices  for  bin  and  supercell  size.  The  main 
focus  of  this  paper  is  on  the  distribution  of  temperature  and  con¬ 
centration  values.  In  this  section,  we  focus  on  the  concentration 
marginal  PDF,  as  its  convergence  is  more  sensitive.  A  common 
model  for  concentration  PDFs  is  the  beta  distribution,  which  we  use 
as  a  surrogate  to  analyze  supercell  and  bin  selections. 

For  a  mass  fraction  x,  the  beta  distribution  has  a  PDF  j 

B(a:,i5)  with  the  normalizing  constant B(a,/5)  = 

This  distribution  is  completely  characterized  by  its  mean  p  and  variance 
a,  related  to  a  and  by  the  formulas  a  =  (p^-p^-p(7^)/(7^  and 
^  =  (p-2p^+p^-(7^+p(7^)/(7^.  In  Fig.  2  left,  we  plot  the  histogram  of  all 
concentration  values  taken  from  the  middle  80%  of  the  mixing  zone  at 
t  =  90,  shortly  after  reshock,  for  Re  =  6  x  10^.  In  this  paper,  we  focus  our 
analysis  on  the  middle  80%  of  the  mixing  zone,  in  order  to  reduce  the 
edge  effects  of  the  mixing  zone.  We  determine  the  parameters  a  =  0.360 
and  ^  =  0.286  from  this  data  and  plot  the  corresponding  beta  distribution 
PDF.  In  Fig.  2  right,  we  compare  the  empirical  (solid  curve)  and  theoretical 
(dotted  curve)  CDFs  of  these  two  distributions.  The  Li  difference,  an 
assessment  of  the  model  error,  is  0.074. 

The  supercell  size  regulates  the  amount  of  stochastic  informa¬ 
tion  captured  in  this  description  and  the  level  of  spatial  resolution 
lost  in  the  process.  In  order  to  assess  mesh  convergence  issues,  we 
choose  supercells  large  enough  that  the  statistical  errors  are  small 
relative  to  the  mesh  errors.  Statistical  sampling  errors  are  typically 
(9(n“^/^),  where  n  is  the  number  of  cells  in  a  supercell. 

We  begin  with  the  model  beta  distribution,  determined  by  the 
parameters  specified  above,  and  take  repeated  random  samples  of 
a  fixed  size.  We  analyze  each  individual  sample,  by  binning  the  data 


Fig.  2.  Left.  Histogram  for  concentration  data  and  beta  distribution  with  the  parameters  a  =  0.360,  ^  =  0.286  (dotted  curve).  Data  from  Re  =  6  x  10^  coarse  grid  (II)  simulation, 
t  =  90.  Right.  Empirical  CDF  of  the  data  (solid  curve)  in  middle  80%  of  the  mixing  zone  from  the  same  simulation  as  left  and  beta  distribution  CDF  (dotted  curve)  with  the  same 
parameters  a, 
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Table  2 

Statistical  convergence  rates  for  repeated  finite  samples  of  fixed  size  drawn  from  the 
beta  distribution,  parameters  as  above.  10  Bins  used  to  construct  the  CDFs. 


Sample  size 

25 

100 

400 

1600 

Mean  of  Li  norm  error  for  CDFs 

0.067 

0.033 

0.017 

0.009 

STD  of  Li  norm  error  for  CDFs 

0.035 

0.018 

0.009 

0.005 

and  forming  a  sample  CDF  from  the  sample  data  and  then  taking  an 
Li  norm  between  the  sample  CDF  and  the  beta  distribution  itself. 
The  Li  norms  are  then  averaged  over  the  repeated  samples  to 
determine  an  expected  mean  statistical  error  that  results  from  the 
supercell  process.  In  Table  2,  we  display  the  means  and  the  STDs  of 
the  statistical  errors  that  arise  for  different  sizing  of  supercells.  We 
can  see  that  upwards  of  400  points  per  supercell  are  necessary  to 
reduce  the  statistical  error  to  a  value  less  than  2%.  This  classifies  the 
tradeoff  between  supercell  sizing  (spatial  resolution)  and  accuracy 
(statistical  error)  of  the  method. 

To  verify  the  beta  distribution  as  a  model  for  the  concentration 
data,  we  also  repeat  the  sampling  process  using  the  empirical  data 
from  the  coarse  grid  (11: 400  x  800)  Re^6  x  10^  simulation  (used  to 
calculate  the  beta  distribution  parameters).  In  the  left  frame  of 
Fig.  3,  we  plot  the  average  statistical  errors  found  when  comparing 
the  beta  distribution  to  finite  samples  drawn  from  it.  In  the  right 
frame,  we  plot  the  average  statistical  errors  found  when  comparing 
empirical  CDFs  of  the  raw  data  from  the  middle  80%  of  the  mixing 
zone  and  random  samples  drawn  from  this  data.  Comparing  the 
differences  between  the  theoretical  formulation  and  our  specific 
application,  we  see  roughly  unchanged  errors  between  the  beta 
distribution  and  the  empirical  distribution,  thus  justifying  the 
choice  of  the  beta  distribution  as  a  surrogate  for  the  concentration 
PDFs.  We  also  show  in  Fig.  3  insensitivity  of  the  statistical  errors  to 
bin  size. 

From  the  above  analysis,  we  draw  the  following  conclusions: 

1.  As  a  result  of  the  insensitivity  of  the  bin  size  parameter,  we  use 
a  bin  sizing  of  10,  per  variable,  in  all  analysis  performed  in  this 
paper. 

2.  To  obtain  statistical  errors  significantly  smaller  than  mesh  er¬ 
rors,  we  choose  to  use  a  supercell  size  close  to  400  data  points 
in  each  supercell  on  the  coarse  grid  (11:  400  x  800). 

3.  High  Reynolds  number  asymptotics 

Flere  we  show,  in  a  purely  hydro  context,  the  mildness  of  the 
variation  introduced  in  passing  from  an  experimental  regime  of 
Re  ~  3.5  X  10"^,  well  studied  experimentally  for  related  Rayleigh— 
Taylor  instabilities,  to  the  ICF  regime  of  Re  ~  6  x  10^  and  larger. 


Table  3 

Li  norm  convergence  of  joint  CDFs  for  temperature  and  con¬ 
centration  (central  column)  and  concentration  only  (right  col¬ 
umn)  for  an  RM  instability  with  Re  =  6  x  10^.  Meshes:  1 
200  X  400,  11  400  X  800,  111  800  x  1600,  IV  1600  x  3200. 


Mesh 

Temp  and  cone 

Cone 

1-lV 

0.092 

0.124 

11-lV 

0.064 

0.079 

111-lV 

0.045 

0.058 

This  fact  supports  our  V&V  proposal  of  experimental  validation  at 
Re  ~  3.5  X  10"^  with  verification  extrapolation  to  Re  ~  6  x  10^. 


3.1.  Convergence  of  joint  temperature-concentration  PDFs 

We  begin  by  analyzing  the  convergence  properties  of  the  joint 
temperature-concentration  CDFs  within  the  middle  80%  of  the 
mixing  zone  for  the  Re  =  6  x  10^  simulation.  Li  norm  convergence 
of  the  joint  temperature-concentration  CDFs  was  established  in 
Ref.  [10]  for  lower  Re  values,  so  this  value  for  Re  represents  an 
extension  of  the  convergence  regime  as  well  as  an  extension  to  the 
case  Pr  =  10“^. 

For  the  supercell  gridding,  since  our  mixing  zone  is  radial,  we 
use  supercells  in  r,  6  space.  Since  the  data  is  homogeneous  in  the  6 
direction  and  sensitive  to  the  r  direction,  discussed  in  further  detail 
in  Sec.  3.2,  we  use  10  supercells  in  the  radial  direction  and  2 
supercells  in  the  angular  direction,  with  a  domain  consisting  only  of 
the  data  points  contained  in  the  common  mixing  zone.  With  the 
data  structured  in  x,y  space,  the  r,  6  bins  have  an  unequal  number  of 
points,  with  an  average  of  642  points  within  each  supercell  of  the 
coarse  grid  (11: 400  x  800).  The  results  of  this  analysis  are  presented 
in  Table  3,  from  which  we  see  evidence  of  norm  convergence  for  the 
CDFs  under  mesh  refinement.  The  slow  convergence  rate  observed 
here  is  consistent  with  the  results  of  [10],  which  used  an  alternate 
convergence  analysis  method. 

We  investigate  the  dependence  of  the  convergence  analysis  on 
the  supercell  griddings  in  Table  4.  We  compare  an  8  x  1  (8  in  r,  1  in 
^),  16  X  2, 32  X  4  and  64  x  8  supercell  grid  structure.  As  the  region 
contained  within  a  supercell  increases  (the  number  of  points 
within  a  supercell  increases,  such  as  moving  from  64  x  8  towards 
8x1),  the  average  norm  error  decreases.  The  average  number  of 
points  in  a  supercell  on  the  coarse  grid  is  approximately  1600, 400, 
100  or  25,  respectively.  We  attribute  the  increase  in  the  norm  error 
to  the  statistical  errors  found  in  Sec.  2.2.  Consistent  with  this 
interpretation  is  the  fact  that  the  increase  of  error  with  supercell 
size  is  nearly  independent  of  the  computational  rows.  We  also  note 
that  the  slight  difference  between  the  1x1  and  1x8  supercell 
gridding  is  due  to  the  convergence  of  the  statistical  errors,  within 


Fig.  3.  Li  norm  of  the  differences  of  the  concentration  CDFs.  Left:  Repeated  finite  samples  of  fixed  size,  drawn  from  the  beta  distribution.  Right:  The  empirical  CDF  from  raw  data  in 
the  middle  80%  of  the  mixing  zone  and  repeated  random  samples  from  this  data.  The  four  sample  sizes  plotted  are:  25, 100,  400  and  1600. 
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Table  4 

Variation  in  the  Li  norm  comparison  of  joint  CDFs  across  mesh  levels  for  different 
choices  of  supercells.  A  coarser  supercell  mesh  corresponds  to  larger  supercells,  and 
improved  statistical  convergence.  All  supercell  choices  are  given  in  radial  x  angular 
direction.  Meshes  as  above. 


Supercell  grid 

1  X  1 

8  X  1 

16  X  2 

32  X  4 

64  X  8 

1-lV 

0.079 

0.085 

0.094 

0.102 

0.138 

11-lV 

0.057 

0.059 

0.065 

0.080 

0.103 

111-lV 

0.039 

0.041 

0.047 

0.059 

0.076 

Fig.  4.  Density  plot  for  Re  =  6  x  10^  after  reshock. 


a  fraction  of  a  percentage  point.  The  supercell  size  (10  x  2)  selected 
is  a  tradeoff  between  statistical  accuracy  and  spatial  resolution. 

3.2.  Spatial  variation  of  mix 

In  Fig.  4,  we  plot  the  density  over  the  mixing  zone  after  reshock, 
for  Re  ^  6  X  10^. 

The  statistics  of  the  PDFs,  as  they  depend  on  x,  y  or  r,  6,  appear 
uniform  in  6  but  the  edge  regions  for  r  differ  from  the  central 


region.  The  unmixed  light  fluid  (in  the  mixing  zone)  is  preferen¬ 
tially  located  near  the  inner  edge  and  the  unmixed  heavy  fluid  is 
similarly  noticeable  near  the  outer  edge  of  the  mixing  zone. 

For  Re  ^  oo ,  we  plot  the  PDFs  in  four  separate  zones  for  the  r 
values,  with  two  central  zones  and  two  edge  zones,  see  Fig.  5.  We 
observe  radial  dependence  in  the  mixing  properties,  with  edge 
zone  PDFs  containing  substantial  amounts  of  unmixed  fluids.  The 
outer  region  has  signiflcant  unmixed  heavy  fluid  and  the  inner 
region  signiflcant  unmixed  light  fluid.  The  central  zones,  especially 
the  outer  central  zone,  are  nearly  equiprobable  in  their  distribution, 
meaning  that  any  mixture  fraction  is  approximately  as  probable  as 
any  other.  The  inner  central  zone  and  the  inner  edge  zone  also  have 
signiflcant  unmixed  light  fluid  concentrations. 

We  conclude  that  the  temperature  and  concentration  fluc¬ 
tuations  are  essentially  stationary  statistically  in  6  and  nearly  so 
in  r  except  near  the  mixing  zone  edges.  We  choose  supercells 
aligned  with  r,  6  coordinate  boundaries,  large  (a)  in  view  of  the 
near  stationarity  of  the  statistics  and  (b)  in  view  of  the  desire  to 
observe  convergence  under  mesh  reflnement,  and  thus  to  reduce 
the  statistical  sampling  errors  relative  to  the  mesh  errors.  As 
with  any  convergence  issue,  the  resolution  will  improve  as  the 
mesh  is  reflned.  Our  intention  is  to  use  stochastic  convergence 
for  assessing  TN  burn.  TN  burn  is  also  a  stochastic  variable,  but 
the  available  observations  are  somewhat  coarse  grained,  and 
appear  to  be  not  greatly  more  reflned  than  the  supercells  in  use 
here. 

3.3.  Re  dependence  of  mix 

As  a  main  result  of  this  section,  we  display  in  Fig.  6,  the  joint 
PDFs  of  temperature  and  concentration  for  the  medium  grid  (111: 
800  X  1600)  Re  =  3.5  x  10"^  and  for  Re  ^  6  x  10^.  Visually  there  is 
minimal  difference  between  the  left  frame,  the  experimental 
regime  with  Reynolds  numbers  of  3.5  x  10^  and  the  right  frame. 


Fig.  5.  PDFs  for  concentration  and  temperature  for  Reynolds  number  Re  =  oo  with  the  variation  of  PDFs  with  the  radius.  Top  row,  left  to  right,  outer  edge  r  values,  and  outer  central 
region.  Bottom  row,  left  to  right,  inner  central  region  for  r  and  inner  edge  region. 
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Fig.  6.  Change  of  Re,  Re  =  3.5  x  10"^  compared  to  i?e  =  6  x  10^.  Left:  Joint  PDF  of  temperature  and  concentration  at  Re  =  3.5  x  10"^.  Right:  Joint  PDF  of  temperature  and  concentration 
at  Re  =  6  X  10^. 


Max:  0.165  f  M  M  ■ 

Min:  0.019  LB  ■  ■ 

Fig.  7.  Li  norm  CDF  differences  between  Fig.  6,  left,  and  right.  The  Li  norms  result  from 
binning  onto  supercells,  to  achieve  numerical  approximation  of  the  associated  CDFs. 

beyond  the  ICF  Reynolds  number  of  6  x  10^.  We  use  a  supercell  grid 
of  10  X  2  as  above,  (approximately  642  points  within  coarse  grid 
supercell),  and  consider  only  the  concentration  marginals.  The  Re 
related  CDF  norm  difference  of  0.077  is  1.3  times  the  mesh  error  of 
0.058  (grids  111-lV)  and  thus  only  marginally  observable.  See  Fig.  7. 


4.  PDF  sensitivity:  an  n  +  1 -parameter  family  of  high  Re 
limits 

The  n-species  compressible  Navier-Stokes  equations,  in  the 
high  Re  limit,  eliminates  Re  as  a  parameter.  There  remain  n  +  1 
dimensionless  fluid  transport  parameters.  The  n  -  1  independent 
Schmidt  numbers,  the  Prandtl  number,  and  the  ratio  of  two  inde¬ 
pendent  viscosities  remain  to  characterize  solutions  of  the  Euler 
equations  achieved  in  the  high  Re  limit  process.  For  LES  it  is 
expected  that  the  limit  Re  ^  oo  is  independent  of  the  choice  of 
physical  transport  parameters,  if  these  are  flxed  and  independent  of 
Re  in  the  inflnite  Re  limit  process.  In  this  sense  the  limit  is  universal 
relative  to  molecular  fluid  transport. 

The  high  Re  limit,  however,  is  sensitive  to  the  path  taken  in  the 
space  of  dimensionless  total  (turbulent  and  molecular)  transport 
coefficients,  if  these  coefficients  are  given  Re  dependent  values.  We 
demonstrate  here  sensitivity  of  the  high  Re  limit  to  the  dimen¬ 
sionless  turbulent  transport  (SGS)  terms  and  their  coefficients.  For 
high  Re  flows,  all  turbulent  Schmidt  numbers  and  the  Prandtl 
number  should  coincide.  There  is  one  parameter  to  set  this  com¬ 
mon  value,  which  is  constrained  by  validation,  i.e.,  comparison  to 
experiment.  A  second  parameter  is  the  ratio  of  bulk  to  shear  tur¬ 
bulent  viscosities.  This  important  point  of  solution  sensitivity  arises 
in  practice  as  dependence  of  the  solutions  upon  the  algorithms  in 
the  high  Re  limit.  The  algorithmic  dependence  of  the  solution  un¬ 
derscores  the  importance  of  validation  and  the  relevance  within 
the  experimental  range  of  Re,  and  a  mild  Re  dependent  and 


theoretically  based  extrapolation,  as  far  as  pure  hydro  issues  are 
concerned,  for  Re  values  beyond  this  range. 

From  this  discussion,  we  see  the  importance  of  turbulence 
models.  A  similar  discussion  (not  included  here)  would  show  the 
importance  of  initial/boundary  conditions  or  forcing  functions.  If 
explicit  SGS  terms  are  not  used,  then  the  flxed  point  is  sensitive  to 
algorithmic  details.  If  a  numerical  algorithm  has  signiflcant 
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Fig.  8.  PDFs  for  the  joint  concentration-temperature  distribution  for  two  values  of  the 
turbulent  transport  coefficient  for  concentration  diffusion.  Above:  10  times  nominal 
and  below:  0.1  times  nominal.  Re  =  6  x  10^. 
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numerical  diffusion  (e.g.,  if  tracking  is  not  used),  then  algorithmic 
issues  will  affect  the  fixed  point  even  if  SGS  terms  are  included. 
Diffusion,  once  added  numerically,  cannot  be  removed  through 
addition  of  SGS  terms.  Thus,  we  anticipate  that  numerical  algorithms 
with  significant  species  concentration  diffusion  will  reach  only 
a  fraction  of  the  full  n  +  1  dimensional  parameter  space  of  high  Re 
limits.  Mesh  refinement  short  of  a  DNS  simulation  will  not  change 
this  statement. 

We  justify  the  assertion  of  solution  nonuniqueness  by  showing 
that  modifications  to  the  dimensionless  turbulent  mass  diffusion 
up  or  down  by  factors  of  10  change  the  PDFs  of  the  solution  vari¬ 
ables.  We  plot  in  Fig.  8,  the  joint  PDFs  for  concentration  and  tem¬ 
perature,  after  an  increase  and  a  decrease  of  turbulent 
concentration  diffusion  by  factors  of  10.  We  see  that  the  high  Re 
solution  is  sensitive  to  the  choice  of  turbulence  model.  Comparing 
spatially  dependent  CDFs  obtained  by  binning  solution  values  into 
supercells,  we  find  an  Li  norm  difference  of  0.187  between  these 
two  simulations,  when  using  a  10  x  2  supercell  grid  and  analyzing 
only  the  concentration  marginal.  This  Li  norm  difference  value  is 
more  than  three  times  as  large  as  the  mesh  error  difference  for  this 
grid  resolution  and  more  than  double  the  normed  difference  be¬ 
tween  the  simulations  with  Re  =  3.5  x  10"^  compared  to 
Re  ^  6  X  10^,  presented  in  Sec.  3.3.  Thus  we  conclude  that  the 
differences  in  turbulent  mass  diffusion  are  significant  and  observ¬ 
able,  in  that  they  yield  CDF  differences  3.2  times  larger  than  the 
statistical  and  mesh  errors. 

5.  Conclusions 

The  results  presented  here  are  a  demonstration,  within  the 
present  research  framework,  that  in  the  high  Re  limit  of  turbulent 
mixing,  atomic  mixing  properties  are  sensitive  to  models  for  tur¬ 
bulent  transport,  and  specifically  to  the  numerical  algorithm.  The 
limit  is  demonstrated  to  be  nonunique,  supporting  the  idea  that 
there  is  an  n  +  1 -parameter  family  of  fixed  points  for  n-species 
mixing,  with  the  fixed  points  labeled  by  the  turbulence  model  or  by 
the  numerical  algorithm.  Among  these  nonunique  choices,  the  FT/ 
LES/SGS  algorithm  selects  a  unique  limit,  which  agrees  well  with 
experiment  within  the  experimental  range  of  Re,  in  contrast  to 
most  ILES  simulations,  which  do  not  [29].  The  turbulent  transport 
coefficients,  within  this  dynamic  SGS  method,  are  uniquely  speci¬ 
fied  by  theory;  the  simulations  have  no  tunable  parameters.  Since 
zero-parameter  numerical  models  of  complex  physical  phenomena 
(i.e.  fluid  turbulence  and  mixing  in  the  present  context)  are  not  so 
common,  we  emphasize  this  aspect  of  the  FT/LES/SGS  algorithm. 
The  algorithm  is,  in  the  opinion  of  the  authors,  unique  among 
compressible  LES  simulations  in  its  level  of  validation  (exper¬ 
imental  confirmation)  for  RT  experiments.  Because  of  the  sensi¬ 
tivity  of  the  atomic  level  mixture  to  turbulence  models  and 
numerical  algorithms,  experimental  validation  is  important. 

A  major  conclusion  of  the  present  paper  is  that  the  ICE  regime  of 
perhaps  Re  ~  6  x  10^  is  only  a  mild  perturbation  (in  terms  of 
strictly  hydro  issues)  of  an  experimental  regime  of  Re  ~  3.5  x  10^, 
for  which  there  are  numerous  hydro  instability  experiments. 

The  main  purpose  of  this  paper  is  to  establish  scientific  results 
which  support  a  two  step  route  with  experiment  and  theory  sup¬ 
ported  extrapolation  to  V&V  for  ICE  simulations. 

The  method  is  suitable  to  the  three  main  instability  stages 
(ablation,  RM  interface  instabilities  and  RT  deceleration)  of  ICE.  We 
explain  the  main  difficulty  (nonuniqueness  of  solutions)  and  its 
resolution  (dynamic  SGS  models  for  LES  and  front  tracking).  We 
review  previous  work  on  the  experimental  comparison  step  and  we 
present  new  material  relating  to  the  extrapolation  step. 

In  future  work,  we  will  examine  this  strategy  in  simulations 
more  closely  modeling  ICE.  In  essence,  our  proposal  is  that  ICE 


hydro  codes  should  be  tested  (validated)  in  the  experimental 
regime  and  tested  numerically  (verified)  for  a  parameter  free 
extrapolation  from  there  to  the  ICE  regime  and  beyond.  Initial 
conditions  for  the  experimental  regime  are  generally  not  known, 
but  we  have  shown  that  these  can  be  reconstructed  from  the  early 
time  data,  with  reconstruction  uncertainty  quantified  and  an 
overall  effect  of  ±5%  on  the  value  of  the  RT  growth  rate  a  [22]. 

This  is  not  the  first  time  in  which  nonuniqueness  has  been  an 
essential  feature  of  the  solutions  of  time  dependent  equations, 
modeled  at  the  inviscid  (Euler  equation)  level.  Shock  refraction 
problems,  which  describe  self  similar  time  dependent  solutions, 
sometimes  have  multiple  solutions.  Elame  speed  is  ambiguous 
when  analyzed  at  the  level  of  the  Euler  equations,  with  the  ambi¬ 
guity  removed  by  consideration  of  the  Prandtl  number.  Detonation 
waves  allow  multiple  solutions,  as  weak  or  strong  detonations  [30], 
while  some  dissipative  mechanism  often  selects  the  weak  deto¬ 
nation.  Equations  suggested  by  three  phase  flow  for  petroleum 
reservoirs  show  nonuniqueness  for  wave  interaction  (Riemann) 
problems  [31  ]. 

Within  the  study  of  turbulence,  the  sensitivity  of  solutions  to 
turbulence  models  is  widely  understood.  Flowever,  for  turbulence 
modeling,  the  analogous  resolution  of  ambiguity  (specification  of 
transport  coefficients)  requires  specifying  the  turbulent  transport, 
exactly  the  quantity  which  introduces  the  ambiguity. 

Flere  we  come  to  the  conundrum  which  is  the  central  theme  of 
this  paper.  Nonuniqueness  of  solutions  is  removed  by  specification 
of  turbulent  (not  laminar)  transport  coefficients,  quantities  not 
accessible  to  direct  measurement  in  a  nearly  infinite  Reynolds 
number  regime,  and  regarding  which  there  is  substantial  dis¬ 
agreement.  The  dynamic  subgrid  models  provide  a  zero  parameter 
solution.  In  other  words  they  do  select  all  needed  coefficients,  and 
remove  all  ambiguity.  The  dynamic  subgrid  models,  combined  with 
front  tracking  provide  excellent  agreement  with  experimental  RT 
data  (validation).  Here  we  claim  that  validation  followed  by  a  mild 
perturbative  verification  step  reaches  the  ICE  regime.  In  view  of  the 
importance  of  this  claim,  we  believe  it  should  be  subject  to  further 
testing  and  studies,  a  task  we  hope  to  carry  out  in  the  future. 
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